
#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include "blaswrap.h"
#include "f2c.h"
#include "clapack.h"
#include "libcell_clapack.h"

void clapack_svd(double* A,int m,int n,double* s,double* uu,double*vt)
{
    char JOBU,JOBVT;
    integer M=m,N=n,LDA=M,LDU=M,LDVT=N,LWORK,INFO;
    integer mn=min(M,N),MN=max(M,N);

    double *wk=(double*)malloc(sizeof(double)*300*MN);
    JOBU='A';JOBVT='A';LWORK=280*MN;
    //LWORK=-1;
    dgesvd_(&JOBU,&JOBVT,&M,&N,A,&LDA,s,uu,&LDU,vt,&LDVT,wk,&LWORK,&INFO);

    free(wk);
}
//*******
//@n: 矩阵A的秩
//返回结果直接储存在A里
//*******
void clapack_inverse(double*A,integer n)
{
    integer *IPIV =(integer*)malloc(sizeof(integer)*n);
    integer LWORK=n*n;
    double* WORK=(double*)malloc(sizeof(double)*LWORK);
    integer INFO;
    dgetrf_(&n,&n,A,&n,IPIV,&INFO);
    dgetri_(&n,A,&n,IPIV,WORK,&LWORK,&INFO);


    free(IPIV);
    free(WORK);
}
//@k:: k是A矩阵的列数

double*  clapack_mult(double* A,double*B, int m,int n,int k)
{
    char transa='N',transb='N';
    integer M=m,N=n,K=k,LDA=M,LDB=K,LDC=M;
    double alpha=1.0,beta=0.0;
    double* re=(double*)malloc(sizeof(double)*M*N);
    dgemm_(&transa,&transb,&M,&N,&N,&alpha,A,&LDA,B,&LDB,&beta,re,&LDC);

    return  re;

}
double* clapack_mult1(double*A,double*B,int m,int n,int k)
{
    char transa='T',transb='T';
    integer M=m,N=n,K=k,LDA=K,LDB=N,LDC=M;
    double alpha=1.0,beta=0.0;
    double* re=(double*)malloc(sizeof(double)*M*N);
    dgemm_(&transa,&transb,&M,&N,&N,&alpha,A,&LDA,B,&LDB,&beta,re,&LDC);


    return  re;
}
